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A key question in evaluation of computer models is Does the 
computer model adequately represent reality! A six-step process for 
computer model validation is set out in Bayarri et al. [Technometrics 
49 (2007) 138-154] (and briefly summarized below), based on com- 
parison of computer model runs with field data of the process being 
modeled. The methodology is particularly suited to treating the ma- 
jor issues associated with the validation process: quantifying multiple 
sources of error and uncertainty in computer models; combining mul- 
tiple sources of information; and being able to adapt to different, but 
related scenarios. 

Two complications that frequently arise in practice are the need to 
deal with highly irregular functional data and the need to acknowl- 
edge and incorporate uncertainty in the inputs. We develop method- 
ology to deal with both complications. A key part of the approach 
utilizes a wavelet representation of the functional data, applies a hier- 
archical version of the scalar validation methodology to the wavelet 
coefficients, and transforms back, to ultimately compare computer 
model output with field output. The generality of the methodology 
is only limited by the capability of a combination of computational 
tools and the appropriateness of decompositions of the sort (wavelets) 
employed here. 
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The methods and analyses we present are illustrated with a test 
bed dynamic stress analysis for a particular engineering system. 

1. Introduction. 

1.1. Validation framework. Bayarri et al. [4] described a general frame- 
work for validation of complex computer models and applied the framework 
to two examples involving scalar data. The framework consists of a six-step 
procedure to treat the combined calibration/validation process, and to as- 
sess the possible systematic differences between model outcomes and test 
outcomes (so-termed biases), by estimating these biases along with uncer- 
tainty bounds for these estimates. The six steps are (1) defining the problem 
(inputs, outputs, initial uncertainties); (2) establishing evaluation criteria; 
(3) designing experiments; (4) approximating computer model output; (5) 
analyzing the combination of field and computer run data; (6) feeding back 
to revise the model, perform additional experiments, and so on. Bayarri et 
al. [3] generalized this work to the situation of smooth functional data arising 
within a hierarchical structure. 

This paper is an extension of Bayarri et al. [4], motivated by methodologi- 
cal needs in analyzing a computer model for analyzing stress on an engineer- 
ing system, a vehicle suspension system, subject to forces over time (the test 
bed problem is described in the following section). This involves the follow- 
ing important — and technically challenging — problems that greatly widen 
the applicability of the six-step strategy. 

Problem 1 (Irregular functional output). In Bayarri et al. [4], exam- 
ples involved real- valued model outputs; in Bayarri et al. [3], outputs were 
smooth functions that could be handled by the simple device of including 
the function argument (time) as an input of the system. In many engineer- 
ing scenarios with functional output, functions are not smooth, and adding 
the function variable to the list of inputs can result in a computationally 
intractable problem. This is so in the test bed problem, for instance, with 
typical irregular functional data indicated in Figure 1. 

Problem 2 (Uncertainty in inputs). A second ubiquitous problem in 
engineering scenarios is that (unmeasured) manufacturing variations are 
present in tested components; incorporating this uncertainty into the anal- 
ysis can be crucial. 

Problem 3 (Prediction in altered or new settings). The point of com- 
puter modeling in engineering contexts is typically to allow use of the com- 
puter model to predict outcomes in altered or new settings for which no field 
data are available. We consider several approaches to this problem. 



COMPUTER MODEL VALIDATION 



3 



1.2. Background and motivation. General discussions of the entire val- 
idation and verification process can be found in Roache [22], Oberkampf 
and Trucano [20], Cafeo and Cavendish [5], Easterling [9], Pilch et al. [21], 
Trucano, Pilch and Oberkampf [26] and Santner, Williams and Notz [24]. 
We focus here on the last stage of the process: assessing the accuracy of the 
computer model in predicting reality, and using both the computer model 
and field data to make predictions, especially in new situations. 

Because a computer model can virtually never be said to be a completely 
accurate representation of the real process being modeled, the relevant ques- 
tion is "Does the model provide predictions that are accurate enough for the 
intended use of the model?" Thus predictions need to come with what were 
called tolerance bounds in Bayarri et al. [4], indicating the magnitude of the 
possible error in prediction. This focus on giving tolerance bounds, rather 
than stating a yes/no answer as to model validity, is important for several 
reasons: (i) Models rarely give highly accurate predictions over the entire 
range of inputs of possible interest, and it is often difficult to character- 
ize regions of accuracy and inaccuracy; (ii) The degree of accuracy that is 
needed can vary from one application of the computer model to another; 
(iii) Tolerance bounds account for model bias, the principal symptom of 
model inadequacy-accuracy of the model cannot simply be represented by 
a variance or standard error. 

The key components of the approach outlined here are the use of Gaussian 
process response-surface approximations to a computer model, following on 
work in Sacks et al. [23], Currin et al. [8], Welch et al. [28] and Morris, 
Mitchell and Ylvisaker [16], and introduction of Bayesian representations 
of model bias and uncertainty, following on work in Kennedy and O'Hagan 
[14] and Kennedy, O'Hagan and Higgins [15]. A related approach to Bayesian 
analysis of computer models is that of Craig et al. [7], Craig et al. [6] and 
Goldstein and Rougier [11, 12], which focus on utilization of linear Bayes 
methodology. 

1.3. Overview. Section 2 describes the test bed example and associated 
data. Section 3 has the formulation of the statistical problem and assump- 
tions that we make for the analysis. In Section 4 we set down the methods 
of analysis with the results in Section 5. Most of the details of computations 
are relegated to the Appendices. 

2. The test bed. The test bed case study is about predicting loads re- 
sulting from stressful events on a vehicle suspension system over time, for 
example, hitting a pothole. In the initial part of the study there are seven 
unmeasured parameters of the system with specified nominal (mean) val- 
ues x nom (referred to later on as Condition A) and unknown manufacturing 
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Table 1 

I/U Map. Uncertainty ranges for the two calibration parameters and the seven input 
parameters that are subject to manufacturing variation 



Parameter Type Uncertainty range 



Damping 1 Calibration [0.125, 0.875] 

Damping^ Calibration [0.125, 0.875] 

xi Nominal + Variation [0.1667,0.8333] 

x 2 Nominal + Variation [0.1667,0.8333] 

x 3 Nominal + Variation [0.2083, 0.7917] 

x 4 Nominal + Variation [0.1923,0.8077] 

x 5 Nominal + Variation [0.3529, 0.6471] 

x 6 Nominal + Variation [0.1471,0.8529] 

x 7 Nominal + Variation [0.1923,0.8077] 



variations S* . There are other relevant parameters that are known and fixed 
and hence not part of the experiments. 

Field data are obtained by driving a vehicle over a proving ground course 
and recording the time history of load at sites on the suspension system. The 
curves must be registered (Appendix A) to assure that peaks and valleys 
occur at the same place. 

In addition, there is a computer model aimed at producing the same re- 
sponse. The computer model is a so-termed ADAMS model, a commercially 
available, widely used finite-element based code that analyzes complex dy- 
namic behavior (e.g., vibration, stress) of mechanical assemblies. The com- 
puter model has within it two calibration parameters u* = (-uj,^) quanti- 
fying two different types of damping (unknown levels of energy dissipation) 
that need to be estimated (or tuned) to produce a matching response. 

For proprietary reasons the specific parameters are not fully described — 
they include characteristics of tires, bushings and bumpers as well as vehicle 
mass. In addition, the values assigned to these parameters are coded on a [0, 
1] scale and the output responses are also coded. In the coded scale the fixed 
values of x nom are all 0.50. The uncertainty ranges for the nine parameters 
were elicited through extensive discussion with engineers and modelers; they 
are given in Table 1, the so-termed Input /Uncertainty (I/U) map (Bayarri 
et al. [4]). Along with the ranges, prior distributions were elicited for (u*, 5*) 
in Section 3.4. 

Field data. In the initial study with Condition A inputs, a single vehicle 
was driven over a proving ground course seven times. The recorded field data 
consist of the time history of load at two sites on the suspension system. 
Plots of the output for Site 1 can be seen in Figure 1 for two of the time 
periods of particular interest. Thus there are seven replicates and a single 
Xnom in the field data. 
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Fig. 1. Model output (bottom) and registered field output (top) for Site 1 at Region 1 
(left) and Region 2 (right). Vertical lines indicate the reference peak locations used in the 
curve registration discussed in Appendix A. 

Computer model runs. A typical model run for the test bed example 
takes one hour, limiting the number of runs that can feasibly be made. 
To select which runs to make we adopted the design strategy used in Ba- 
yarri et al. [4]: the 9-dimensional rectangle defined by the ranges of the 
parameters in Table 1 is first transformed into the 9-dimensional unit cube 
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and a 65-point Latin hypercube design (LHD) is then selected using code 
by W. Welch that finds an approximately maximin LHD. The actual de- 
sign points, as well as a 2-dimensional projection of the design, can be 
found in Bayarri et al. [2]. In addition, we added a point at the center 
(0.5, . . . , 0.5), the nominal values. One run failed to converge and was deleted 
from the experiment leaving a total of 65 design points. 

3. Formulation, statistical model and assumptions. 

3.1. Formulation. Given a vector of inputs, x = (x±, . . . ,Xd), to a time- 
dependent system, denote the "real" response over time t as y R (x;t). Field 
measurement of the real response has error and we write the rth replicate 
field measurement as 



where the e r (-)' s are independent mean zero Gaussian processes (described 
below). Some inputs may also have error; we take that into account below. 

In addition, there is a computer model aimed at producing the same 
response. The computer model may have within it calibration/tuning pa- 
rameters u = (ui, . . . ,u m ) that need to be estimated (or tuned) to produce 
a matching response. The model output is then of the form y M (x,u;t); it 
is affected by u but the real response, y R , is not. The connection between 
model output and reality is then expressed in 



where u* is the true value of the (vector) calibration parameter, y (x, u*;i) 
is the model response at time t and the true value of u, and 6(x;i), defined 
by subtraction, is the associated bias. In situations where u is a tuning 
parameter, there is no "true value" so u* should be thought of as some type 
of fitted value of u, with the bias defined relative to it. 

Data from the field and from the computer model runs provide informa- 
tion for estimating the unknowns in (1) and (2). The Bayesian analysis we 
employ takes note of the fact, as in Bayarri et al. [4], that the unknowns u* 
and the bias are not statistically identifiable and, consequently, specification 
of their prior distributions is of particular importance for the analysis. 

Inputs. Some inputs may be specified or physically measured with es- 
sentially perfect accuracy. Those that remain fixed for both field data and 
model runs play no further role and are not part of x. Other (unmeasured) 
inputs will have specified nominal values (generally, they will vary in the 
experiments) that are subject to manufacturing variation with specified dis- 
tributions. We write these as 

(3) x = x nom + d, 






y H (x;t)=y M (x,u*;t) + b(x;t) 
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where x nom is the known nominal value and the distribution of the manu- 
facturing variation 8 can be specified. In effect this transforms (2) into 

y R (x nom + S*;t) = y M (x nom + S*,u*;t) + 6(x nom + 6*;t), 

where S* is the actual (unknown) value of S. The parameters 3 are like 
calibration parameters in that they are unknown but physically real. 

The unknowns. Prior to making computer runs or collecting field data, 
the unknowns in (1) and (2) are (y , u*,S*, b, V e ), where V e is the covariance 
function of e. A full Bayesian analysis would contemplate placing priors on 
these unknowns and, given field data and model runs, produce posterior 
distributions. But the complexity (e.g., of irregular functional output) and 
high dimensionality militate against such a strategy unless simplifications 
can be made. One such is the use of a basis representation of the functional 
data. In particular, to handle the irregular functions, we will consider wavelet 
decompositions. Other settings may allow different representations such as 
Fourier series or principal components (Higdon et al. [13]). 

3.2. Wavelet decomposition. The nature of the functions in Figure 1, 
for example, suggests that wavelet decomposition would be a suitable basis 
representation (see Vidakovic [27]; Miiller and Vidakovic [18] and Morris et 
al. [17] are among other references with applications related to ours). 

The wavelet decomposition (more details are in Appendix A) we use for 
y M is of the form 

(4) y A/ (x,u;t)=^ U ; l M (x,u)^(t), 

i 

where the wavelet basis functions ^i(t) are default choices in R wavethresh 
(Daubechies wavelets of index 2; for simplicity of notation, we include the 
scaling function as one of the basis elements). Similarly, the field curves (rth 
replicate) are represented as 

(5) yf(x;t)=$>£(x)*i(*)- 

i 

A thresholding procedure, used to produce a manageable number of coeffi- 
cients while maintaining adequate accuracy, leads to the approximations 

y M (x,u;t)=^ U ;f(x,u)^(t), 

(6) 

The accuracy of the approximations using the reduced set of elements for 
the test bed problem is indicated in Figure 2. Since similarly accurate re- 
construction resulted for all field and model run curves, we also assume that 
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Fig. 2. The original curve and wavelet reconstructed curve for the first field-run at Site 1. 



reality and the bias function can be accurately represented by the same basis 
elements, and write 



y R (x;t) = ^ U ;f(x)^(t) ! 



(7) 



&(x;t) = 5>?(x)*i(f), 



iei 



iei 



Matching coefficients and using (1) and (2), we get 
(8) ™f(x)=™f(x,u*) + ^(x) Vie/, 

(9) 



w( r (x) = Wi t (x) + e 



Vie/. 



That the measurement errors in the wavelet domain, S %y . £1X6 normally dis- 
tributed with mean zero and are independent across replications r follows 
from the Gaussian process assumption on the process measurement error. 
We will, in addition, assume that these coefficients are independent across 
i, with possibly differing variances a\. This independence assumption is im- 
posed so that later computations are feasible. That might seem unrealistic, 
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given that the seven residual functions {yf (x;t) — y^f(x;i)) — shown in the 
top of Figure 3 — can be seen to be correlated in time t, suggesting that the 
ei r should be correlated in i. But, as long as the of differ, independent nor- 
mal ei r lead to a Gaussian process with mean zero and covariance function 
J2iei °1 The bottom of Figure 3 gives seven realizations of this 

process, with the of being estimated by the usual unbiased estimates, based 
on the replicates. The correlation patterns of the two processes seem quite 
similar. 

Our approach is to analyze each of the retained wavelet coefficients in (8) 
and (9), and recombine them to obtain estimates and uncertainties for the 
"original" functions in (1) and (2). 

3.3. GASP approximation. For y M , the wavelet coefficients are func- 
tions of (x, u). Because we cannot freely run the computer model for every 
(x,u), we approximate each of the retained coefficients using data from 
computer runs. Formally, we start with a Gaussian process prior distribu- 
tion on a coefficient wf^(x, u). Given computer model runs, y M (xfc,Ufc), 
where {(xfe,Ufe),A; = 1,...,K} are the design points in a computer experi- 
ment, we extract the data {wf 1 (x^, u^)} and approximate wf 4 ^, u) as the 
Bayes predictor, the posterior mean, of wf' 1 (x, u) given the data. 

The Gaussian process priors we use are as in the GASP methodology 
described in Bayarri et al. [4]. Let z = (x, u). For each i £ I (the set of 
retained wavelet coefficients), the GASP assumption is that wf 1 (z) is a 
Gaussian process with mean /j, M , constant variance 1/Xf 1 and correlation 
function 

Ml l\ I qM\ i \2-a M \ 

4 (z,z)=exp^2^P%\zp-z p \ » j, 

where um is the number of coordinates in z, the j3's are nonnegative param- 
eters and the a's are between and 1. 

Let Of = {/if, \f,af£,/3fg;p = 1, . . . ,n M } be the collection of the (hy- 
per) parameters determining the Gaussian prior distribution of wf' 1 . To pro- 
duce the Bayes predictor we have to deal with the Of'^s; we do so in Section 
4. 

3.4. Other prior specifications. Priors for u*,<5* are context specific. En- 
gineering advice led to adopting uniform priors for u* on their ranges in 
Table 1. For the manufacturing variations of the unmeasured parameters in 
Table 1, the advice led to normal priors with standard deviations equal to 
1/6 of the ranges of the uncertainty intervals in Table 1. Specifically, 

7r(iii) = ir(u 2 ) = Uniform on [0.125,0.875], 
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Fig. 3. The seven residual field processes (top) and seven simulated error processes (bot- 
tom). 



7r(<5i) = 7r(<5 2 ) ~ iV(0,0.1111 2 ) truncated to [-0.3333,0.3333], 
7r(<5 3 ) ~ iV(0,0.09723 2 ) truncated to [-0.2917,0.2917], 
7r(<5 4 ) = tt(5 7 ) ~iV(0,0.1026 2 ) truncated to [-0.3077,0.3077], 
7r(<5 5 ) ~iV(0,0.04903 2 ) truncated to [-0.1471,0.1471], 
7r(<5 6 ) ~ iV(0,0.1176 2 ) truncated to [-0.3529,0.3529]. 
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The of are given the usual noninformative priors 7r(<r?) oc 1/crf. Priors 
for the bias wavelet coefficients, the (u^)'s, will be Gaussian but restricted 
to depend only on those coordinates of x nom that vary in the experiments. 
Because there is only one x nom in the test bed experiment, we only consider 
the case where the (w^O's are constants, though a similar approach can be 
taken for more general settings. 

In the wavelet decomposition, each i E I belongs to some resolution level, 
(in the test bed — see Appendix A — the levels go from to 12). It is 
natural and common to model wavelet parameters hierarchically, according 
to their resolution level. The priors for the biases w\ at resolution level j are 

(10) 7rK b |r|)~AT(0,rJ), j = 0, . . . , 12. 

This is a strong shrinkage prior, shrinking the biases to zero. One might be 
concerned with such strong shrinkage to zero, but the computer modeling 
world is one in which biases are typically assumed to be zero, so that utilizing 
a strong shrinkage prior has the appeal that any detected bias is more likely 
to be believed to be real in the community than would bias detected with a 
weaker assumption. (Of course, there are also statistical arguments for using 
such shrinkage priors.) 

The hypervariances rj are assigned a variant of a typical objective prior 
for hypervariances, 



j 

where crj = average of of for i at level j. The <r| provide a "scale" for the 
Tj and are necessary — or at least some constants in the denominators are 
necessary — to yield a proper posterior. 

4. Estimation and analysis. We restrict attention for the most part to 
the context of the test bed. It will be clear that much of what is done can 
be generalized. 

Approximating wavelet coefficients. In Appendix A we find that, for the 
test bed, there are 289 wavelet coefficients wf 4 to be treated. The Gaussian 
prior for each has 20 hyperparameters (coordinates of Of 1 ). A full Bayesian 
treatment would then require treatment of 5780 parameters, an infeasible 
task. Instead, we treat each coefficient separately and estimate each Of 1 by 
maximum likelihood based on the model-run data = {wf 1 (z^)}, using 
code developed by W. Welch. Recall that z = (x, u) and denote the fcth 
design point in the computer experiment by Zfc. For the test bed there are 
65 Zfc's. 

Letting i = {jlf 4 , Xf 4 , af£ , Mf ; p = 1,...,um} be the maximum likeli- 
hood estimates of the O s, it follows that the GASP predictive distribution 
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of wf 4 (z) at a new z is 

(11) wf (z) I wf,df ~N(m? (z),^ M (z)), 

where 

^ M (z) = ft" + 7f (z)'(f f - /if 1), 
^ M (z) = ^ 7 -7f( Z ) / (ff)" 1 7 i M (z), 

with 1 the vector of ones, I\ (65 x 65 in the test bed example) the co- 
variance matrix for the model-run data w" estimated by plugging in Q i , 
and 

7 J M (z) = (l/Af)(cf / (z 1 ,z),...,cf(z fc ,z))', 

where c is the estimated correlation function. For the rest of the paper we 
will use (11) as the definition of the GASP predictive distribution. 

A plot of the (2,/3) pairs can be found in Bayarri et al. [2]. A J3 near zero 
corresponds to a correlation near one and a function that is quite flat in that 
variable; an a near zero corresponds to a power of two in the exponent of 
the correlation, suggesting smooth functional dependence on that variable. 
For most of the pairs (5, (3), one or the other is near zero. 

Full justification of the use of the plug-in maximum likelihood estimates 
for the parameters 6^ is an open theoretical issue. Intuitively, one ex- 
pects modest variations in parameters to have little effect on the predictors 
because they are interpolators. In practice, "Studentized" cross-validation 
residuals (leave-one-out predictions of the data normalized by standard er- 
ror) have been successfully used to gauge the "legitimacy" of such usage (for 
examples and additional references see Schonlau and Welch [25] and Aslett 
et al. [1]). Only recently, Nagy, Loeppky and Welch [19] have reported sim- 
ulations indicating reasonably close prediction accuracy of the plug-in MLE 
predictions to Bayes (Jeffreys priors) predictions in dimensions 1-10 when 
the number of computer runs = 7 x dimension. 

The posterior distributions. Restricting to the test bed problem, we sim- 
plify the notation by referring only to S, the deviation of x from the nominal 
inputs x nom , and rewrite (8) and (9), Mi £ I, as 

12 

where the independent N(0,af). 

The field data can be summarized by the (independent) sufficient statistics 

w[ = i £ W ( r (6*), sj = £«r(<n - *Bf y. 

r=l r=l 
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(We drop the argument d* for these statistics, since the statistics are actual 
numbers given from the data.) Key facts to be retained, using (11) and (12) 
and properties of normal distributions are 



u 



2 2 



(13) wf | wf (6*,u*)M, °t ~ N(< (**, u*) + < K)' 

where (11) is used to get the last expression. 

Let w 6 , r 2 , <r 2 , d* and u* denote the vectors of the indicated parameters 
and write their prior distribution as 

/ b 2 2c* * \ 

tt(w ,T ,(T ,6 ,u ) 

: 7T(W 6 I T 2 ) X 7r(<5* , U* , T 2 | <T 2 ) X 7r(cr 2 ) 
7 2 12 



(14) 



n-K 6 ir? (l) ) 



i=o 



li=l 



i=l 



iei 

The data, from field and computer model runs, can be summarized as D = 
{ujf,s?,w^;i = l,...,289}. Using (11), (13), (10) and (14), together with 
standard computations involving normal distributions, it is straightforward 
to obtain the posterior distribution of all unknowns as 

W?A^U*),W^*,U*,<T 2 ,T 2 | D) 

= W^ A/ (^>*)!w 6 ,<r,u*,<x 2 ,D) 

(15) X7r post (w 6 |5*,u*,o- 2 ,r 2 ,D) 

X TT post (S\u*,T 2 | cr 2 ,D) 
X TTpost(o- 2 | D), 



where 



7r post (w M (S*, u*) | w fc , 6*, u*, <x 2 , D) - J] JV(m u , ^); 



mit 



^ M (<T,u*) 



(16) 



^(<P,u*) + (l/7K? 
(V7)a 2 



+ 



^(r,u*) + (i/7)^ 



(mf(<T,u*)) 
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It 



^(<5*,u*)(l/7)a! 



VV(S*,u*) + (l/r)a. 



and 



(17) 



vr post (w 6 | <T,u*,<t 2 ,t 2 ,D) ~Y[N(mx,V2i); 

-(<-mf(5*,u*)), 



m2i 



V: 



2i 



^(r,u*) + (l/7K? + r| w 

^(^(J'.iQ + qA)*?) 

F.M (<5 * )U * ) + (1/7K 2 + ^ ( . ) - 



The third factor in (15) is 

TT post {S*,\l*,T 2 | <T 2 ,D) 

oc L(W F ,s 2 I <T,u*,cr 2 ,-r 2 ) x 7r(<T,u*,T 2 | a 2 ), 



(18) 



where the marginal likelihood L found by integrating out w 6 and w (8*, u*) 



in the product of the full likelihood and 7r(w | 



L(w J 



c* „* _2 _^ 
, U , <T , T 



n 



' 6 1 • 
i 



x cxp 

Finally, the fourth factor in (15) is 

1 



L e ;^(r,u*) + (i/7K 2 +4) 

1/ (<-mf(<T,u*)) 2 



7rp OS t{(T I D) OC 



(19) 



n 



iei 
x 



RF exp 



2V^(r,u*) + (l/7)a 2 +r 2 w 



2a 



|L(w F ,- 



s 2 |<r,u*,<r 2 ,T 2 )d«5*du* dr 2 . 



At this point we make an approximation, and ignore the integral in (19); 
that is, we simply utilize the replicate observations to determine the poste- 
riors for the a 2 . The reason for this is not computational; indeed, one can 
include the af in the posterior in (18) and deal with them by a Metropolis 
algorithm. Instead, the motivation is what we call modularization, which is 
meant to indicate that it can be better to separately analyze pieces of the 
problem than to perform one global Bayesian analysis. The difficulty here 
is that there is a significant confounding in the posterior distribution be- 
tween the calibration parameters, the bias function, and the af and this, for 
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instance, allows bias to be replaced by larger af. Here we have seven repli- 
cate observations for each of, so simply utilizing the replicate observation 
posteriors and preventing the confounding has intuitive appeal. (Formaliz- 
ing this argument is not so easy. The difficulty occurs because part of the 
model — the modeling of the bias — is quite uncertain. Better or more robust 
modeling of the bias may correct the problem within a full Bayesian analy- 
sis, but the difficulty of doing so argues for the simpler modular approach. 
We will discuss these issues more fully elsewhere.) 

Simulating from (16) and (17) and the first factor of (19) is, of course, 
trivial, but simulating from (18) requires MCMC methodology. Given the 
complexity of the problem, the MCMC requires careful choice of proposal 
distributions in order to achieve suitable mixing. Discussion of these pro- 
posals is relegated to Appendix B because of the level of detail needed to 
describe them, but we note that these are technically crucial for the method- 
ology to work and required extensive exploration. 

The end result of the simulation is a sample of draws from the posterior 
distribution in (15): each saved draw from the first factor of (19) is used to 
generate the MCMC sample for the third factor, with both being used to 
generate a draw from the second and the first factors, using (16) and (17). 
We saved every 200th draw from 200,000 MCMC iterations for the third 
factor, thereby obtaining a final sample of 1000 draws, 

(20) {w M ' h (S* h , u* h ),w bh ,d* h , u* h , a 2h , r 2h ; h = l,..., 1000}. 
The results below are based on this sample from the posterior. 

5. Results. 

5.1. Estimates of S*,u*. Histograms for d*,u* (Figure 4) are obtained 
by forming a histogram for each component of S* , u* from the corresponding 
elements in (20). The calibration parameters are moderately affected by the 
data but, of the input variables, only x§ and xq have posteriors that are 
significantly different from the priors. The posterior for X5 is piled up at the 
end of the allowed range for the variable, which suggests the (undesirable) 
possibility that this uncertain input is being used as a tuning parameter to 
better fit the model; a case could be made for preventing this by additional 
modularization. 

5.2. Estimation of bias and reality. Posterior distributions of the bias 
and reality curves are obtained by recombining the wavelets with the pos- 
terior wavelet coefficients from (20). For instance, the posterior distribution 
of b is represented by the sample curves 

(21) b h {t) = Y,w\ h %(t), /i = l,...,1000. 
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. 4. Histogram of the posterior draws for the two calibration parameters (ui and ui ), 
the seven input parameters (xi through xr), along with their priors (solid lines). 
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Fig. 5. The estimate of the bias function (solid line) with 90% tolerance bounds (dashed 
lines) for suspension Site 1 and at Region 1. 

The posterior mean curve, b(t) = X)Ji=i b h (t), is plotted as the solid 
line in Figure 5. The uncertainty of this estimate of b is quantified by pro- 
ducing upper and lower uncertainty (tolerance) bounds at each t by, for 
example, taking the lower a/2 and upper 1 — a/2 quantiles of the posterior 
distribution of b(t), that is, 

L h (t) = - quantile of {b h (t); h = l,..., 1000}, 

(22) 

U b (t) = (l-^j quantile of {b h {t); h = l,..., 1000}. 

These lower and upper bounds are also plotted in Figure 5. It is apparent 
in Figure 5 that the bias function is significantly different from 0, especially 
in the neighborhood of 8.7 and 9.1. 

The bounds in (22) are symmetrically defined. Alternative tolerance bounds 
can be defined by only requiring that 100a% of the curves lie outside the 
bounds; a useful choice would satisfy this condition and minimize the height 
of the bounds, U b (t)-L b (t). 

Figures 4 and 5 provide marginal distributions of u* , x nom + S* and the 
bias, but it is important to note that these are highly dependent in the 
posterior. Hence most analyses involving these quantities must be based on 
their joint, rather than marginal, distributions. 
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Estimating reality with uncertainty bounds is done similarly: take the 
sample of wavelet coefficients w Rh = wf /Ih (S* h , u* h ) + wf 1 and form 

y Rh (t) = Y^^ h Mt), 

i 

1000 

h 

| quantile of {y Rh (t); h = l,..., 1000}, 
(\ - ^ quantile of {y Rh (t); h = l,..., 1000}. 

We call y R (t) the bias- corrected prediction of reality. Figure 6 exhibits the 
bias-corrected prediction and associated uncertainty band. 

Figure 6 further shows a comparison between bias-corrected prediction 
and pure model prediction, the latter being defined as 

(24) f(t)=^mf(U)$,(f), 

i 

where d = j^jY,h^* h an d Q = j^qq J2h u * h an ^ mf 1 {S,^) is the posterior 
mean of the wavelet coefficients with plugged-in estimates for the unknown 
parameters [use (11)]. 

In practice, it may be that running the computer model after estimating 
S* , u* is feasible. Then an alternative (and preferred) pure model prediction 
is y M (S,u;t). 

Assessing the uncertainty for predicting reality by the pure model pre- 
diction [see (24)] can be done by considering samples {y Rh (t) — y (i)} 
and forming bounds. If a new model run producing y M (d,u;t) is used in- 
stead, then the computation should be redone with this additional model 
run included, but a typically accurate approximation is to simply consider 
{y Rh (t) — y M (S,u;t)} and proceed as before. Here it may be useful to con- 
sider asymmetric bounds because the pure model predictions may lie entirely 
above or below the realizations of reality. But plots like that of Figure 6 al- 
ready show the gap between pure model prediction and reality. 

5.3. Predicting a new run; same system, same vehicle (new run). In 
some prediction settings, not necessarily the one of the test bed, there is 
interest in predicting a new field run with the same inputs (and the same 
system). Prediction is done by adding in draws from a N(0,af h ) distri- 
bution and then following the same prescription as in (23) to form y F (t) 
and corresponding uncertainty bounds. This, of course, produces wider un- 
certainty bounds. 



(23) 



y R (t) 

L R (t) 
U R {t) 
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Fig. 6. Bias-corrected prediction of reality (solid black line) with 90% tolerance bands 
(dashed black lines), pure model prediction (solid grey line) and field runs (light grey lines). 



5.4. Extrapolating. There are many follow-on settings where prediction 
is called for. We single out three such: (i) the same vehicle except that 
x nom changes; (ii) same vehicle type but new components, that is, the same 
nominal values and prior distributions for the manufacturing variations but 
a new realization of the <5's from the prior distribution; (hi) a new vehicle, 
different nominal (x-y, . . . ,xj) inputs, but the same prior distributions for 
the <Ts. 

Any analysis we perform for the test bed is severely constrained by the 
fact that we have limited field data on just one set of inputs measured with 
error. 

5.4.1. Same vehicle, different load. Here the same vehicle was tested but 
with added mass. This causes a known change in x nom with everything else 
(including the <5's) remaining unchanged. The effect of a simple change in 
inputs such as this can be addressed using the difference in computer model 
runs at the nominal values of the inputs. More formally, suppose we add a 
computer model run at (x nom + A, u nom ) where A is the modest change in 
the nominal inputs (A has only one nonzero coordinate if the only change is 
in the mass). Suppose we also have a run at the old nominals (x nom , u nom ); 
if not, use the GASP prediction based on the existing runs. Our assumption 
then is that (with x* referring to the true unknown input values for the 
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original system) 

y M (x* + A,u*;t)-y M (x*,u*;t) 

— V (^nom ~i~ A, U nom , t) y (x nom , U nom , t) = D(t). 

We can then make predictions under the new inputs by simply adding D(t) 
to the old predictions. 

This is illustrated in Figure 7; the given bias corrected prediction and 
tolerance bands for the vehicle with the additional mass are simply the 
results of Section 5.3 translated by D(t). The grey solid line is the pure 
model prediction. The light grey line is the actual result from a field test of 
the system with added mass, and the strategy appears to have successfully 
extrapolated in the critical Region 1. Again, the pure model prediction is 
not successful. 

5.4.2. New vehicle of the same type. In this setting the nominal values 
x nom remain the same but the new 5's are random draws from their prior 
(population) distribution and are, therefore, different from those for the 
field-tested system. This is of particular interest in practice, in that it is 
prediction for the population of vehicles of the given type that is of prime 
interest, rather than just prediction for the single system/vehicle tested. 

The calibration parameters u* do not change; they belong to the model 
and, if physically real, are inherently the same for all systems/ vehicles of the 
same type. Denote the parameters of the new components by z new = (x nom + 
<5 new , u* ) . The input parameters of the computer runs remain z^ = (x nom + 
<5fc,Ufc); and z* = (x nom + <5*,u*) are the true values for the tested system. 
Denote the associated model wavelet coefficients for the new components by 
w AI (z ncw ). 

Since 5 ncw is independent of (w b , 5*, u*, a 2 , r 2 , D), the predictive (pos- 
terior) distribution is (with the Lj's denoting likelihood terms arising from 
the data) 

TT pos t(w M (z new ),w M (z*), w fc ,£ ncw ,<T,u*,cr 2 ,-r 2 | D) 
(25) oc4<5 ncw ) x^(w 6 ,<T,u*,cr 2 ,T 2 ) 

x L 1 (w M (z new ),w M (z*),{wf(z k )} | <5 ncw ,<T,u*) 
xL 2 ({^},{ S 2 }|u; M (z*),w b , < T 2 ) 
oc7r poS i(u; A/ (z ncw ) I w M (z*),{wf (z k )},5 ncw ,S* ,u*) 
x 7r(<5 ncw ) x 7r(w b ,<T,u*,cr 2 ,i- 2 ) 
xL 3 (u; M (z*),K f (z fc )}|^,u*) 
xi 2 ({™,},{ Sl 2 }|/(z'),wV 2 ). 
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To sample from (25), note that the last three factors in the expression 
yield exactly the same posterior for (w M (z*),w b ,<5* ,u*,er 2 ,r 2 ) as before 
(with the same modularization used), so the draws from the previous MCMC 
can be used in the new computations. Since S new can be drawn from its prior 
7r(<$ ncw ), it only remains to draw from TT post {w M (z) | w M (z*), {wf 1 (z k )}, S new , 
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(5*,u*). But this is simply the GASP distribution where w (z*) has to be 
added to the model run data. Therefore, one simply determines the GASP 
for the augmented runs wf = [wf 1 {z\ ) , wf 1 (z2 ) , • • • , wf 1 (z*. ) , wf 1 (z* ) ) , that 
is, 

(26) (sw) | wf°, Of ~ iV(mf /0 (z ncw ), t/ 4 M0 (z ncw )), 

-A/ 

with #j as in Section 4 and 

mf °(z new ) = /2f + 7 f °(z new )'(f f V(wf - /If 1), 

where 7f°(zncw) = (1/Af )(cf (z 1? z new ), . . . , cf (z fe , z ncw ), cf(z*, z ncw ))' and 
Y i is obtained by appending the column 7f °(z*) and the row 7f °(z*)' 

tor, . 

Application of these expressions yields draws h = 1, . . . , 1000 from the 
posterior distribution of the ith wavelet coefficient for the new system as 

w[ h {z h ) = wf(z h ) + wf + 4, 

where e\ ~ iV(0, of 1 ). Fi gure 8 plots the predictions of the new system with 
uncertainty bands. The uncertainty has increased as compared with Figure 
6 not only because the prior for <5 new is used rather than the posterior for the 
tested vehicle, but also because these are predictions of field runs (instead 
of reality). 

5.4.3. New vehicle with new nominals. The primary engineering use of 
computer models is to extrapolate to a system with new nominals when 
there are no new field data. This will require strong assumptions, especially 
about the bias. The simplest assumption about the bias, which we make 
here, is that the new system has the same bias function as the old system. 
The calibration parameters u* are also assumed to remain the same. We use 
the joint — and highly dependent — posterior distribution of the bias and u* 
from the original system extensively in what follows. 

The new system has the same I/U map as the original system, but with 
new nominal values x„ om (which we refer to as Condition B). The new £'s are 
taken to have the same priors as for Condition A. The same 65-point design 
on (<5,u) was used as before with the addition of one central value. Again 
one run failed, leaving 65 usable model runs. Registration of the output 
curves is unnecessary because there are no field data and the computer runs 
are assumed to be inherently registered. The resulting computer runs were 
passed through the same wavelet decomposition as before, retaining only 
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Fig. 8. Predictions for a new vehicle of the same type. The solid black line is bias-cor- 
rected prediction, the dashed black lines are 90% tolerance bands, the grey line is pure 
model prediction and the light grey lines are field runs. 



those coefficients that appeared earlier. The resulting GASP for the ith 
wavelet coefficient is 



W: 



.BM , 



■) 



W; 



.BM 



N(fhf M (z),Vr VI (z)), 



BM, 



(27) 

exactly as in (11). 

This new GASP analysis is done only with Condition B data. A GASP 
analysis combining the model runs for Condition B and the model runs for 
the original system is not used because the changes in the nominals are too 
large to safely assume connections between the computer model runs for the 
two systems. 

The situation now is analogous to that of the previous argument for new 
vehicles of the same type with the same nominals. In the current case, again 
using the independence of S B from the other unknowns, the predictive (pos- 
terior) distribution of the relevant unknowns can be written as 



7T(w BM (z B ),W fe ,^,U 



a ,t 

b S B 



D) 



7r(w BA/ (z B ) | w°, 5^,11 , cr- 
xtt{S b ) xtt(w 6 ,uV 2 ,t 2 I 

vr(w BA/ (z B )|<5 B ,u*,D) 
xtt{S b ) X7r(w 6 ,u*,cr 2 ,T 2 



r 2 ,D) 



D) 
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here 7r(w BA/ (z B ) | <5 B ,u*,D) is just the GASP distribution in (27) and 
tt(S b ) is the prior for B inputs from the I/U map. Draws of w 6 ,u*,<r 2 ,r 2 
are made from the old posterior distribution for the original system. Because 
w b and u* are highly dependent in the posterior, they must be jointly sam- 
pled for the extrapolation; naive approaches-such as simply trying to add 
the bias from Figure 5 to the pure model prediction-will not succeed. 

The "carry-over" assumptions for the bias and the field variances lead to 
draws from the posterior distribution of the wavelet coefficients for B to be 



w BFh {z B h)=w BM {z h )+w bk + £ 



where 4 ~ N(0, of 1 ). 

In the top of Figure 9 the prediction for B is presented. Actual field data 
(8 replicate runs for B) were available afterward (not used in constructing 
the predictions and tolerance bands) and they are superimposed (in light 
grey) on the plots in the top of Figure 9; darker grey curves represent model 
runs. The effectiveness of carrying over the assumptions from A to B is 
apparent. 

If a strong assumption, such as the constancy of bias, is to be made it is 
best to be extremely careful about implementing the assumption. Here, for 
instance, physics considerations might suggest that an assumption of con- 
stant multiplicative bias is more sensible than an assumption of constant 
additive bias. (Use of a multiplicative bias is mentioned in Fuentes, Guttorp 
and Challenor [10].) A standard way of implementing multiplicative bias is 
to repeat the above analysis using the log of the various curves. A simpler 
alternative is to transform the additive biases obtained above into multi- 
plicative biases, and apply these multiplicative biases to the GASP draws 
under Condition B. Bias in the additive system can be written 

b h (t) = y Rh (t)-y Mh (t); 

the corresponding multiplicative representation of the bias is 

b h (t) yRk{t) 1 

"multW — yMhU\ ' 

which would lead to draws from the posterior for reality under Condition B 
of 

y BR h{t) = y BM h{t)x{1 + h h^ m 

The bottom panel of Figure 9 presents the analogue of the top using the 
multiplicative bias. The additive and multiplicative predictions are not much 
different. The next section discusses another site in the suspension system 
for which the analysis in the paper was implemented, called Site 2. For this 
site, the multiplicative predictions under Condition B were noticeably better 
than the additive predictions, as indicated in Figure 10. 
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Fig. 9. Prediction at Site 1 of a new vehicle under Condition B in Region 1. Top: 
additive bias. Bottom: multiplicative bias. The solid black line is bias-corrected prediction, 
the dashed black lines are 90% tolerance bands and the grey lines are model runs; the light 
grey lines are the field runs later provided and not used in the analysis. 

6. Site 2 analyses. Analyses for Site 2 of the system proceed in exactly 
the same way as those for Site 1. The posterior distributions for the cal- 
ibration parameters (u\,U2) as well as for 8 are in Figure 11. These are 
somewhat different than those for Site 1 in Figure 4. These parameters are, 
of course, the same for either site, but the limited data available at each 
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FlG. 10. Prediction at Site 2 of a new vehicle under Condition B in Region 1. Top: 
additive bias. Bottom: multiplicative bias. The solid black line is bias- corrected prediction, 
the dashed black lines are 90% tolerance bands and the grey lines are model runs; the light 
grey lines are the field runs later provided and not used in the analysis. 

site lead to somewhat different posterior distributions. In particular, X5 no 
longer seems to be used as a tuning parameter, but it is possible that 112 is 
being so used. A natural solution to help reduce such potential overtuning 
is to do a bivariate functional analysis of the two sites jointly. This is being 
pursued separately. 
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Fig. 12. Bias estimate with 90% tolerance bounds for Site 2 at Region 1. 

Figure 12 presents the estimated bias function and Figure 13 the bias cor- 
rected and pure model predictions, along with the corresponding tolerance 
bounds. The presence of bias can be clearly seen, as well as the noticeably 
improved performance of the bias corrected prediction over the pure model 
prediction. Other figures for Site 2 are omitted since, with the exception of 
Figure 10, they do not provide further insight. 

APPENDIX A: DATA REGISTRATION AND WAVELET 

DECOMPOSITION 

For the wavelet representations of the output curves it is important that 
the same wavelet basis elements simultaneously represent the important 
features of all curves. In the test bed problem the heights of the peaks and 
valleys of the curves from the field data are of primary importance, but their 
locations are not the same across the curves, due to random fluctuations in 
the tests. Thus we first align the curves so that the major peaks and the 
major valleys occur at the same locations. In other applications, alignment 
would likely be based on other key features of the curves. (In the test bed, the 
timing of the major events is not of concern — only the forces at these events 
are of interest. If it were important for the computer model to accurately 
reflect timing, as in the analysis of airbag deployment in Bayarri et al. [3], 
this mode of registration could not be used.) 




We did not try to align the curves from the computer model runs, since 
variation in these curves could not be ascribed to random fluctuation. (One 
might think that the computer model curves would be automatically aligned 
but, surprisingly, in the test bed they did show some misalignment, perhaps 
due to differences in the damping parameters.) We construct a reference 
curve (for alignment of the field curves) by averaging the model-run curves 
and use piecewise linear transformation to align the peaks and valleys of the 
field curves to this reference curve. The details follows: 

Step 1. Construct a dyadic grid (points of the form j/2 q on the interval 
where the function is defined; for the test bed, the interval [0, 65] 
covered the range of importance and q = 12). For each computer 
run extract the output values on the dyadic grid. Construct a pseu- 
dooutput for points not in the grid by linear interpolation and treat 
them as actual outputs. 

Step 2. From the K computer runs {K = 65 in the test bed) define the 
reference curve as y M (t) = Y$=i U M ( x k, Ufc! t). 

• For the first major event, located in the region 6 < t < 11, define 

- A = location (value of t) of the maximum of the reference curve 

y M (t); 

- = location of the maximum of y^; 
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- a = location of minimum of y M (t); 

— af = location of minimum of yf. 

• For the second major event, located in the region 37 < t < 41, 
define B, B^, b, b^ analogously. Assume a < A < b < B. 

Step 3. For each r, match ,Af with a, A by transforming t in [a^, Af] to 

t =a + {t-a r )—p — -p. 

Now define the registered yf on the interval [a, A] as registered y^{t') 
original y^(t), where t' and t are connected as above. 
Step 4. Assume that b < B. As in Step 3, register yf on the interval [Af , bf\ 
by mapping [A^ , bf] into [A, b] via t' = A + (t — A^) & f!^f ■ Similar 

registrations are done for the intervals [0, a^] , [b^ , Bf\ and [B^, 65] . 

Figure 1 shows the registered data for Site 1 of the suspension system. 
"Time" in the figures is not literal time, but is a convenient scaling of the 
realigned time. 

The wavelet decompositions are described in part in Section 3.2. The ba- 
sis functions are such that at "level" there is a scaling constant and for 
level j = 1, . . . , 12 there are 2- ?_1 basis functions. To balance the need for 
approximation accuracy with the need to minimize the number of terms for 
computational feasibility, we considered each model-run and field curve and 
retained all coefficients at levels through 3; for levels j > 3, we retained 
those coefficients that, in magnitude, were among the upper 2.5% of all co- 
efficients at all levels for the given function, according to the R wavethresh 
thresholding procedure. We then took the union of all resulting basis func- 
tions for all the model-run and field curves. For the test bed there were 231 
retained elements for the output from Site 1 on the system, and 213 for 
the output from Site 2. The combined (from both sites) number of retained 
elements was 289 and we used these for all analyses. The indices attached 
to these 289 retained basis elements are denoted by /. 

APPENDIX B: THE MCMC ALGORITHM 
Step 1. For h = 1, . . . , 1000, sample the of 1 from the distribution 

/ 2 \ / 2 
InverseGamma ^3, —^J , ^shape = 3, scale = 

Step 2. For h = l, . . . , 1000, make draws S* h , u* h , T 2h from the posterior dis- 
tribution in (18). (This is complicated — the process is described 
last.) 

Step 3. Given 6* h ,u* h ,a 2h ,r 2h , draw w bh from the distribution in (17). 
(This is simply done by making a draw, for each i, from a normal 
distribution with the specified means and variances.) 
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Step 4. Given S* h ,u* h ,cr 2h ,T 2h ,w bh , make a draw of w Mh from the distri- 
bution in (16). (Again this is simply done by draws from normal 
distributions.) 

For Step 2, we use a Metropolis-Hastings scheme to generate the (h + l)st 
sample. We break this up into two substeps. 

Step 2.1. Propose r 2 by generating from q(r 2 \ r 2h ) = Yli= qi(T 2 \ r 2h ), 
where 

, { — if T 2 G \ T 2h e -0.7 2h 0.71 

(28) ft (7f |rf)oc r 2 ' T 6 ,T ' 6 J ' 

[ 0, otherwise. 

The posterior density of r 2 is not very spiked, so this type of fairly 
broad local proposal works well. Finally, form the Metropolis- 
Hastings ratio 

_ ir(S h ,u h ,r 2 \a 2h ,-D)q(T 2h \t 2 ) 
P ~ 7r(S h ,u h , r 2h | a 2h , D) q(r 2 \ T 2h ) 

and define r 2 ( h+1 ) = r 2 with probability min(l,p); r 2 ^ +1 ) = r 2h 
otherwise. 

Step 2.2. Let T% = [a 5 k ,A 5 k ] and T% = [a£, AJJ] denote the intervals on which 
the prior densities for the corresponding variables are nonzero, and 
define 

n Sh = [max(4,^ -0.05), min(4,4 l + 0.05)], 
T^ uh = [max(aj|,ii£ - 0.05), min(A£,u£ + 0.05)]. 
Propose S, u from 

g(d,u\d\u h ) 

= {{{\U{5 k \Ti) + \U{5 k \Tf h )) 

k=l 

*Y\{\u{u k \m + \u(u k \Tz uh )). 

k=l 

The logic here is that the posterior densities for some of the param- 
eters are quite flat, so that sampling uniformly over their support 
(Tfc or Tfc) would be quite reasonable as a proposal. On the other 
hand, some of the posteriors are quite concentrated, and for these 
it is effective to use a locally uniform proposal, centered around 
the previous value and with a maximum step of 0.05; this leads 
to uniforms on Tf h or T£ uh , which are the regions defined by the 
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intersection of the local uniforms and the support of the priors. 
Since the goal here was to create a procedure that can be auto- 
matically applied for this type of problem, 50-50 mixtures of the 
two proposals were adopted. 

Finally, form the Metropolis-Hastings ratio 



and set (d^ h+1 \ u( h+1 )) = (S,u) with probability min(l,p), and 
equal to (S h ,u h ) otherwise. 

These Metropolis-Hastings steps typically yield highly correlated iterations, 
so we actually cycle through them 200 times (with fixed <r 2h ) before saving 
the variable values for feeding into Steps 3 and 4. 
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